#include "dtcl4.hpp"

dtcl4::dtcl4(matrix<double> par){
  params=par;
}

dtcl4::dtcl4(const dtcl4& d){
  if(this != &d){
    params=d.params;
  }
}
  
dtcl4& dtcl4::operator=(const dtcl4& d){
  if(this != &d){
    params=d.params;
  }
  return *this;
}


vector<double> dtcl4::operator()(double t){
  
 unsigned n=params.size1();
  vector<double>r(n);
  for(unsigned ii=0;ii<n;ii++){
    double delta,gamma,lambda;
    gamma=params(ii,0);
    lambda=params(ii,1);
    delta=params(ii,2);
    double C=t*lambda;
    double A=(gamma*lambda)/(pow(lambda,2)+pow(delta,2));
    double B=pow(gamma,2)*pow(lambda,5)*exp(-C)/(2*pow((pow(lambda,2)
							+pow(delta,2)),3));
    double D=delta*t;
    double E=delta/lambda;
    
    r(ii)=A*(1-exp(-C)*(cos(D)-E*sin(D)))
      +B*((1-3*pow(E,2))*(exp(C)-exp(-C)*cos(2*D))
	  -2*(1-pow(E,4))*C*cos(D)
	  +4*(1+pow(E,2))*D*sin(D)
	  +E*(3-pow(E,2))*exp(-C)*sin(2*D));
  }
  return r;
}

